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Abstract. 

We study the topological susceptibility, x, in QCD with two quark flavours using lattice 



<L> 
> 

field configurations that have been produced with an 0(a) improved quark action. We 
find clear evidence for the expected suppression at small quark mass, and examine the 
variation of \ with this mass. The resulting estimate of the pion decay constant, /„. = 
105 ± 6 IJo MeV, is consistent with the experimental value of ~ 93 MeV. We compare 
X to the large-iV c prediction and find consistency over a large range of quark masses. 
We discuss the benefits of the non-perturbative action improvement scheme and of the 
stategy of keeping the lattice spacing (nearly) fixed as the quark mass is varied. We 
compare our results with other studies and suggest why such a quark mass dependence 
has not always been seen. 



1 Introduction 



In gluodynamics (the pure gauge or "quenched" theory) lattice calculations of the con- 
tinuum topological susceptibility now appear to be relatively free of the systematic errors 
arising from the discretisation, the finite volumes and the various measurement algorithms 
employed (for a recent review, see U). 

The inclusion of sea quarks in ( "dynamical" ) lattice simulations, even at the relatively 
large quark masses currently employed, is numerically extremely expensive, and can only 
be done for lattices with relatively few sites (typically 16 3 32). To avoid significant finite 
volume errors, the lattice must then be relatively coarse, with, in our case, a spacing 
a ~ 0.1 fm. This is a significant fraction of the mean instanton radius, as calculated 
in gluodynamics, and thus precludes a robust, detailed study of the local topological 
features of the vacuum in the presence of sea quarks. The topological susceptibility, on 
the other hand, may be calculated with some confidence and provides one of the first 
opportunities to test some of the more interesting predictions of QCD. Indeed, it is in 
these measurements that we find some of the most striking evidence for the presence of 
the sea quarks (or, alternatively, for a strong quenching effect) in the lattice simulations. 

We recall that the ensembles used here have been produced with two notable features 
[0. The first is the use of an improved action, such that leading order lattice discreti- 
sation effects are expected to depend quadratically, rather than linearly, on the lattice 
spacing (just as in gluodynamics). In addition, the action parameters have been chosen 
to maintain a relatively constant lattice spacing, particularly for the larger values of the 
quark mass. 

These features have allowed us to see the first clear evidence for the expected sup- 
pression of the topological susceptibility in the chiral limit, despite our relatively large 
quark masses. From this behaviour we can directly estimate the pion decay constant 
without needing to know the lattice operator renormalisation factors that arise in more 
conventional calculations. 

The structure of this paper is as follows. In Section |2| we discuss the measurement 
of the topological susceptibility and its expected behaviour both near the chiral limit, 
and in the limit of a large number of colours, N c . In Section [| we describe the UKQCD 
ensembles and the lattice measurements of the topological susceptibility over a range of 
sea quark masses. We fit these with various ansatze motivated by the previous section. 
We compare our findings with other recent studies in Section |j. Finally, we provide a 
summary in Section [5|. 

These results were presented at the IOP2000 @], the Confinement IV gj and, in a much 
more preliminary form, the Lattice '99 || conferences. Since then, we have increased the 
size of several ensembles and included a new parameter set to try and address the issue of 
discretisation effects in our results. We also have more accurate results from the quenched 
theory with which to compare. A brief summary of this work appears in [Q. 
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2 The topological susceptibility 



In four-dimensional Euclidean space-time, SU(3) gauge field configurations can be sep- 
arated into topological classes, and moving between different classes is not possible by 
a smooth deformation of the fields. The classes are characterised by an integer-valued 
winding number. This Pontryagin index, or topological charge Q, can be obtained by 
integrating the local topological charge density 

q(x) = \e^ T F^{x)F a aT {x) (1) 

over all space-time 

Q = ^ J q(x)d i x . (2) 

The topological susceptibility is the expectation value of the squared charge, normalised 
by the volume 

X=—- (3) 

An isolated topological charge induces an exact zero-mode in the quark Dirac operator. 
As a result sea quarks in the vacuum induce an instanton-anti-instanton attraction which 
becomes stronger as the quark masses, m u , m^, . . . , decrease towards zero (the 'chiral 
limit'), and the topological charge and susceptibility will be suppressed to leading order 
in the quark mass 0, 

x=s(— +— y 1 (4) 

\m u m d J 

where 

S = - lim lim (Ol^lO) (5) 

Vtlq— »Q V^oo 

is the chiral condensate (see |7| for a recent discussion). Here we have assumed (Ol^VlO) = 
(0|mt|0) = (0\dd\0) and we neglect the contribution of heavier quarks. PCAC theory 
relates this to the pion decay constant f n and mass M n via the Gell-Mann-Oakes-Renner 
relation as 

f*M> = (m u + m d )Z + 0(m 2 q ) (6) 
and we may combine these for Nf degenerate light flavours to obtain 

f 2 M 2 

X = ^ + 0{Mt) (7) 

in a convention where the experimental value of f n ~ 93 MeVQ. This relation should hold 
in the limit 

flMlV^l, (8) 

1 N.B. there is a common alternative convention, used for internal consistency in some more preliminary 
presentations of this data 0], where a factor of 2 is removed from in Eqn. [?], and where f„ is a factor 
of a/2 larger, around 132 MeV. 
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which is the Leutwyler-Smilga parameter to leading order in M%. We anticipate our 
results here to say that even on our most chiral lattices the lhs of Eqn. ^|is of order 10 and 
so this bound is well satisfied. Thus a calculation of x as a function of M v should allow us 
to obtain a value of /„-. This method has an advantage over more conventional calculations 
in that it does not require us to know lattice operator renormalisation constants which 
are required for matching matrix elements, but which are usually difficult to calculate. In 
principle, we require instead knowledge of the renormalisation of the topological charge 
operators, but we see in the next Section that this problem can be readily overcome. 

As m q and M n increase away from zero we expect higher order terms to check the rate 
of increase of the topological susceptibility so that, as m q , M n — > oo, x approaches the 
quenched value, x qu - I n fact, as we shall see below, the values of x that we obtain are not 
very much smaller than x qu . So there is the danger of a substantial systematic error in 
simply applying Eqn. [7| at our smallest values of M n in order to estimate /„■. To estimate 
this error it would be useful to have some understanding of how x behaves over the whole 
range of m q . This is the question to which we now turn. 

There are two quite different reasons why x might not be much smaller than x qn - The 
obvious first possibility is that m q is large. The second possibility is more subtle: m q may 
be small but QCD may be close to its large- N c limit ||. Because fermion effects are non- 
leading in powers of N c , we expect x ~^ X qu f° r any fixed, non-zero value of m q , however 
small, as the number of colours N c — > oo. There are phenomenological reasons || [| for 
believing that QCD is 'close' to N c = oo, and so this is not an unrealistic consideration. 
Moreover in the case of D = 2 + 1 gauge theories it has been shown |K| that even SU(2) 



is close to SU(oo). Recent calculations in four dimensions [JTTJ, |TJ] indicate that the same 
is true there. In the present simulations, the lighter quark masses straddle the strange 
quark mass and so it is not obvious if we should regard them as being large or small. We 
shall therefore take seriously both the possibilities discussed above. 

We start by assuming the quark mass is small but that we are close to the large-iV c 
limit. In this limit, the topological susceptibility is known to vary as 



n 

where x°° i /oo are the quantities at leading order in N c . In the chiral limit, at fixed 
N c , this reproduces Eqn. [7|. In the large-iV c limit, at fixed M ff , it tends to the quenched 
susceptibility x°° because /£, oc N c . The corrections to Eqn. || are of higher order in M% 
and/or lower order in N c . 

We now consider the alternative possibility: that m q is not small, that higher-order 
corrections to x wm t> e important for most of the values of m q at which we perform 
calculations, and that we therefore need an expression for x that interpolates between 
m q = and m q = oo. Clearly one cannot hope to derive such an expression from first 
principles, so we will simply choose one that we can plausibly argue is approximately 
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correct. The form we choose is 



^J^-pf*'"^) (10) 

where f n is the pion decay constant in the chiral limit. The coefficients have been chosen 
so that this reproduces Eqn. |7] when M n — > and x ~^ X qu + when M n — > oo. 

Thus, this interpolation formula possesses the correct limits and it approaches those limits 
with power-like corrections. 

We shall use the expressions in Eqns. [7], || and [LC] to analyse the m q dependence of 
our calculated values of x an d to obtain a value of f„ together with an estimate of the 
systematic error on that value. In addition, the comparison with Eqn. |9] can provide us 
with some evidence for whether QCD is close to its large- iV c limit or not. 

3 Lattice measurements 

We have calculated x on fi ye complete ensembles of dynamical configurations produced by 
the UKQCD collaboration, as well as one which is still in progress 0. Details of these data 
sets are given in Table [I]. The SU(3) gauge fields are governed by the Wilson plaquette 
action, with "clover" improved Wilson fermions. The improvement is non-perturbative, 
with c sw chosen to render the leading order discretisation errors quadratic (rather than 
linear) in the lattice spacing, a. 

The theory has two coupling constants. In pure gluodynamics the gauge coupling, f3, 
controls the lattice spacing, with larger values reducing a as we move towards the critical 
value at (3 = oo. In simulations with dynamical fermions it has the same role for a 
fixed fermion coupling, k. The latter controls the quark mass, with k —>■ n c from below 
corresponding to the massless limit. In dynamical simulations, however, the fermion 
coupling also affects the lattice spacing, which will become larger as n is reduced (and 
hence m q increased) at fixed (3. 

The three least chiral UKQCD ensembles (by which we mean largest m^/mp) are e6, es 
and e 4 . By appropriately decreasing (3 as n is increased, the couplings are 'matched' to 
maintain a constant lattice spacing |L3|, [14| (which is 'equivalent' to (3 ~ 5.93 in gluody- 
namics with a Wilson action 0) whilst approaching the chiral limit. The physical volume 
and discretisation effects should thus be very similar on these lattices. The remaining en- 
sembles have lower quark masses, but are at a slightly reduced lattice spacing. (To have 
maintained a matched lattice spacing here would have required reducing (3 to values where 
the non-perturbative value c sw (/3) is not known.) As the lattices are all L > 1.5 fm, we 
believe that the minor reduction in the lattice volume should not lead to significant finite 
volume corrections. We also remark that ensembles e2 and e3 have been matched to have 
approximately the same chirality, but at (mildly) different lattice spacings. 

Four-dimensional lattice theories are scale free, and the dimensionless lattice quantities 
must be cast in physical units through the use of a known scale. For this work, we use the 
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Sommer scale |L5| both to define the lattice spacing for the matching procedure, and to 



set the scale. The measured value of fo on each ensemble, as listed in Table [I], corresponds 
to the same physical value of r = 0.49 fm. (f is the dimensionless lattice value of r in 
lattice units i.e. f = r (a)/a. We use the same notation for other quantities.) As we are 
in the scaling window of the theory, we can then use the naive dimensions of the various 
operators to relate lattice and physical quantities, e.g. r^x — r o 4 X + C( a2 ); where we 
have incorporated the expected non-perturbative removal of the corrections linear in the 
lattice spacing. 

Further details of the parameters and the scale determination are given in @. Measure- 
ments were made on ensembles of 400-800 configurations of size L 3 T = 16 3 32, separated 
by ten hybrid Monte Carlo trajectories. Correlations in the data were managed through 
jack-knife binning of the data, using ten bins whose size is large enough that neighbouring 
bin averages may be regarded as uncorrelated. 

We begin, however, with a discussion of lattice operators and results in the quenched 
theory. 

3.1 Lattice operators and % qu 

The simplest lattice topological charge density operator is 

<?W = ^^Tr U^(n)U aT (x) (11) 

where £/ M „(n) denotes the product of SU(3) link variables around a given plaquette. We 
use a reflection-symmetrised version and form 

Q = ^X>», d2) 

X = Jjf (13) 

with a 4 (L 3 T) the lattice volume. In general, Q will not give an integer-valued topological 
charge due to finite lattice spacing effects. There are at least three sources of these. First 
is the breaking of scale invariance by the lattice which leads to the smallest instantons 
having a suppressed action (at least with the Wilson action) and a topological charge less 
than unity (at least with the operator in Eqn. |TT|). We do not address this problem in 
this study, although attempts can be made to correct for it fL6] , but simply accept this 
as part of the overall 0(a 2 ) error. In addition to this, the underlying topological signal 
on the lattice is distorted by the presence of large amounts of UV noise on the scale of 
the lattice spacing [O, and by a multiplicative renormalisation factor |T3 that is unity in 



the continuum, but otherwise suppresses the observed charge. Various solutions to these 
problems exist [[J. In this study we opt for the 'cooling' approach. Cooling explicitly 
erases the ultraviolet fluctuations so that the perturbative lattice renormalisation factors 
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for the topological charge and susceptibility are driven to their trivial continuum values, 
leaving 0(a 2 ) corrections that may be absorbed into all the other lattice corrections of 
this order. We cool by moving through the lattice in a 'staggered' fashion, cooling each 
link by minimising the Wilson gauge action applied to each of the three Cabibbo-Marinari 
SU(2) subgroups in the link element in turn. (The Wilson gauge action is the most local, 
and thus particularly efficient at removing short distance fluctuations whilst preserving 
the long range correlations in the fields.) Carrying out this procedure once on every link 
constitutes a cooling sweep (or 'cool'). The violation of the instanton scale invariance on 
the lattice, with a Wilson action, is such that an isolated instanton cooled in this way 
will slowly shrink, and will eventually disappear when its core size is of the order of a 
lattice spacing, leading to a corresponding jump in the topological charge. Such events 
can, of course, be detected by monitoring Q as a function of the number of cooling sweeps, 
n c . Instanton-anti-instanton pairs may also annihilate, but this has no net effect on Q. 
However, these observations do motivate us to perform the minimum number of cools 
necessary to obtain an estimate of Q that is stable with further increasing n c (subject to 
the above). 

To estimate this point we calculate the normalised correlation function between the 
topological charges measured after n c cooling sweeps, and a nominally asymptotic 25 
cooling sweeps: 

(Q(n c )Q(25)) 

Rn[n c ) = — -t—^ 7. r-. (14) 

1 f (Q(n c )) 2 + (Q{25))* 
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In Fig. |l| we show a typical plot for ensemble e 4 . As discussed before, we have opted 
not to attempt to round the topological charge to integer values. We find (Q)(n c ) and 
(Q 2 )(n c ) to be stable within statistical errors for n c > 5, and the results presented here 
are for n c = 10. 

The topological charge of a configuration is related to the smallest eigenvalues of the 
Dirac matrix and as such is often believed to be one of the slowest modes to decorrelate 
during Monte Carlo simulations. It is crucial for the error analysis that the bin sizes for 
the data are at least twice the integrated autocorrelation times. In Fig. ^| we plot a typical 
time series of the topological charge measured every ten hybrid Monte Carlo trajectories. 
The rapid variation between configurations suggests that the integrated autocorrelation 
time is small even for the topological charge. Estimates of this are given in Table ^ in 
units of ten trajectories. These reinforce the impression gained from the time series plots. 
The bins used in the jack-knife statistical analysis are between 400 and 1000 trajectories 
in length and thus may be confidently assumed to be statistically independent. It is 
interesting that although autocorrelation times are hard to estimate accurately, it does 
appear that they increase as we move away from the chiral limit (c.f. Ref. 0). 

In Fig. |3] we divide the topological charge measurements made over an ensemble into 
bins of unit width centred on the integers, and plot a histogram of the probability of 
finding a configuration with each charge, with errors from the jack-knife analysis. We 
find for all our ensembles that these histograms are very close to being symmetric, centred 
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around Q = and consistent with a Gaussian envelope. The hybrid Monte Carlo appears 
to be sampling the topological sectors correctly, and it is legitimate to extract an estimate 
of the topological susceptibility. On the histograms we show this estimate as a Gaussian 
curve 



Q- 



2i>n < 15 > 

The central line uses our estimate of x, whilst the outlying curves use the central value 
plus or minus one standard deviation. The agreement with the histograms is good. 

We also remark that on a lattice one obviously loses instantons with sizes p < 0(a). 
Since the (pure gauge) instanton density decreases as p 6 when p — > this would appear 
to induce a negligible 0(a 7 ) error in the susceptibility. However this is only true for 
a — > 0, and the error can be substantial for the coarse lattices often used in dynamical 
simulations. 

In general, then, we expect the topological charge and susceptibility to be suppressed 
at non-zero lattice spacing. In gluodynamics with the Wilson action this suppression 
can, typically, be fitted by a leading order, and negative, 0(a 2 ) correction term starting 
from quite moderate values of (3. Of course different ways of calculating the topological 
charge differ substantially at finite values of a, even if they agree in the continuum limit. 
(See, for example, Table 27 in |T2|].) An important factor for non-zero a is whether the 
topological charge is rounded to the nearest integer after cooling or not, as can be seen 



in Table 8 in 1 11]. A gluodynamic calculation of x that uses a method very similar to 



the one used in the present paper, in particular an unrounded topological charge, can be 



found in the last column of Table 8 in JTTJ] . The values of the susceptibility listed there 
can be fitted, for (3 > 5.7, by 

r*x {qn) = 0.065 (3) - 0.28 (4)/rg. (16) 

In obtaining this fit we have used the interpolation formula for f ((3) that is given in JT9 



We shall use this formula in Section | as a guide to the typical variation of x with a. 

In the next Section we shall want to compare our calculations of x, as obtained in the 
presence of sea quarks, with an appropriate quenched limit. The 'equivalent' quenched 
limit will, of course, depend on the lattice spacing, i.e. on the value of f . However, 
because of our strategy of varying (3 so as to approximately match the values of f at 
the different values of m q , the variation in this quenched value, f^x^, over the range of 
lattice spacings of our ensembles, is in fact much less than the statistical errors on the 
measurements themselves. (If the calculations had been performed at fixed (3, the lattice 
spacing would have become increasingly coarse with increasing m q , and the reduction in 
the quenched susceptibility would have been much more pronounced over this range of k. 
We shall return to this important point when we discuss other work in section f|.) One 
finds [[| that fo = 4.714 (13) at (3 — 5.93 in the quenched theory, demonstrating that this 
provides an appropriate quenched limit for our calculations (see Table p. At (3 = 5.93 
the interpolation formula for r (j3) that we used to obtain Eqn. [16] gives f — 4.735 
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and if we insert this in Eqn. [16] we obtain r^x^ = 0.0525 (11). Different methods for 
calculating fo differ by 0(a 2 ) terms, and part of the difference between the UKQCD 
value, f = 4.714(13), and the Sommer value, f — 4.735, might be due to this. If we 
rescale the susceptibility to account for this, then the value of r^x^ drops to ~ 0.0516. 
To do better than this we need to take into account the fact that the calculations of the 
topological charge that enter into Eqn. [16] are obtained by methods that are not exactly the 
same as those used in the calculations of the present paper. The potentially significant 
differences are that 20 cooling sweeps and an unsymmetrised topological charge were 



used in ||1 1|| . while we use 10 cooling sweeps and a symmetrised charge. To estimate the 
systematic shift induced by these differences we have performed calculations on 300 16 4 
lattice field configurations generated at (3 = 5.93 (separated by 50 Monte Carlo sweeps). 
We find that there is no significant difference between the susceptibility as calculated 
by the symmetrised and unsymmetrised charges, whether after 10 or 20 cooling sweeps. 
There is, on the other hand, a small but significant difference between the susceptibility 
as calculated after 10 and 20 cooling sweeps. This reduces our estimate of the equivalent 
quenched susceptibility by ~ 0.0025(7). So taking all this into account we take our 
equivalent quenched susceptibility to be given by 

r* X (ciu) (P = 5.93) = 0.049 (2). (17) 



3.2 Sea quark effects in the topological susceptibility 

In Table |^ we give our estimates of the topological susceptibility in physical units, using 
r as the scale. In Fig. |] we plot TqX versus a similarly scaled pseudoscalar meson mass 
(calculated, of course, with valence quarks that are degenerate with those in the sea, 
i.e. ^valence = K sea)- We also plot the corresponding value of the quenched topological 
susceptibility, as calculated at (3 — 5.93. 

Comparing the dynamical and quenched values, the effects of the sea quarks are clear. 
Whilst the measurement on and es are consistent with the quenched value, moving to 
smaller m q (oc Mj) the topological susceptibility is increasingly suppressed. 

We can make this observation more quantitative by attempting to fit our values of f^x 
with the expected functional form in Eqn. [?|, so extracting a value of f n . But we must first 
be clear whether this fit is justified, and what exactly we are extrapolating in, bearing in 
mind that Eqn. is strictly a chiral expansion that describes the behaviour for small sea 
quark masses in the continuum limit. 

An immediate concern is that our cooling technique will occasionally misidentify the 
value of Q and, in addition, that at finite lattice spacing the exact zero modes associated 
with the topological charge Q are shifted away from zero. All this implies that x wm n °t in 
fact vanish as m q — > 0. However we expect this effect to be small for the following reasons. 
First, the use of an improved fermion action should ensure that the zero-mode shift will 
only be significant for very small instantons, i.e. those whose sizes are p ~ 0(a). These 
are unlikely to survive the cooling and should not contribute to our calculated value of Q. 
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Large instantons, on the other hand, for which any zero mode shift should be negligible, 
will certainly survive the cooling. The remaining ambiguity involves the smaller, but not 
very small, instantons. These might be erased by the cooling but the probability is small 



simply because the number of such charges is small [|16|, ^0|. For example, we can see from 



Fig. 12 in |TTJ that the cooling only appears to cut out instantons with p < 2.5a. Since 
the lattice spacing in that plot is a factor of 1.5 smaller than at our equivalent quenched 
(3 value of 5.93, we would expect that cooling at (3 — 5.93 would affect instantons that 
have p < 4a in that plot. We see from the SU(3) curve in that figure (after scaling up 
by a factor of ~ 4 to take us from a volume of 20 4 at (5 = 6.2 to a volume of 16 3 32 at 
(3 = 5.93) that this involves less than one topological charge per field configuration. This 
is a small effect in the present context. Thus it is reasonable to assume that this effect 
can be neglected for values of a and m q comparable to the ones that we study, and that 
we can then apply Eqn. |7| to values of the susceptibility obtained by varying m q at a 
fixed value of the lattice spacing: except that now the decay constant /„- will be the one 
appropriate to that lattice spacing. 

Now, whilst most of our data points are evaluated on a trajectory of constant lattice 



spacing in the parameter space ||13|| , not all are. If fof^ varied significantly with a over this 
range of a, it would not be clear how to perform a consistent chiral extrapolation through 
the data points. The non-perturbative improvement of the action, however, removes the 
leading order lattice spacing dependence and the residual corrections in this range of lattice 
spacings appear to be small, at least in measurements of (quenched) hadron spectroscopy 



21|) p2| . An indication of the possible size of the effect on topological observables comes 
from comparing our two measurements of f^x & t {roM n ) 2 ~ 4. The range of lattice 
spacings here (f varies from 4.75 to 5.14) is comparable to that over our total data set. 
The accompanying shift in the topological susceptiblity is, however, within our statistical 
errors. Accordingly, we proceed now to attempt a common chiral extrapolation to the 
data, assuming throughout that lattice corrections to the relations discussed before are 
too small to be discernable in our limited data. We shall at the end of this Section return 
to the issue of scaling violations. 



For this purpose, and in the light of the discussion at the end of section pTT] , it is useful 
to redisplay the data in Fig. |5], where the leading order chiral behaviour would then be a 
horizontal line, 

"M = c (18) 
Ml 

and including the first correction gives a generic straight line 



2 -v 

c + Cl (r M n ) 2 . (19) 



Ml 



In each case the intercept is related to the decay constant by c = (roA-) 2 /4. We now 
follow a standard fitting procedure, first using the most chiral points, then systematically 
adding the less chiral points until the fit becomes unacceptably bad. The larger the 
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number of points one can add in this way, the more evidence one has for the fitted form 
and the more confident one is that the systematic errors, associated with the neglected 
higher order corrections, are small. The results of performing such fits are shown in 
Table |3] and those using the two and four most chiral points respectively are plotted on 
Fig. [5]. We see from the Table that the fits using Eqn. [19] show much greater stability and 



these are the ones that will provide our eventual best estimate for f w . 

We should comment briefly on the determination of the fitting parameter errors. In 
performing all but the constant fit we must contend with the data having (small) errors 
on the abscissa in addition to the ordinate. In order to estimate their affect on the fitting 
parameters, we first perform fits to the data assuming that the abscissa data take their 
central values. Identical fits are then made using the central values plus one, and then 
minus one standard deviation. The spread of the fit parameters obtained provides what 
is probably a crude over-estimate of this error (given there is some correlation between 
the ordinal and abscissal uncertainties) but is sufficient to show that it is minor. We 
show this spread as a second error, and for estimates of the decay constant we add it in 
quadrature to the other fit parameter error. 

It is remarkable that we can obtain stable fits to most of our data using just the first 
correction term in Eqn. [19[ Nonetheless, as we can see in Fig. [|, our values of the 
susceptibility are not very much smaller than the M n = 00 quenched value and we need 
to have some estimate of the possible systematic errors that may arise from neglecting the 
higher order corrections that will eventually check the rise in \. As discussed earlier we 
shall do so by exploring two possibilities. One is that the reason why \ is close to x qu is 
not that m q is 'large' but rather that N c = 3 is large. Then the values of \ should follow 
the form in Eqn. [|. A second possibility is simply that our values of m q are indeed large. 
In that case we have argued that the functional form Eqn. [H^ should be a reasonable 
representation of the true mass dependence. We now perform both types of fit in turn. 

We begin with the first possibility, and therefore fit the data with the following ansatz 



j c c 3 (f M 7r ) 



4 = — y , (20) 

c 3 + co(r il4) 2 ' 

where we expect C3 = r 4 x qu up to 0(1/N 2 ) corrections^. To test this we fit up to seven 
data points. The first six are measured in the dynamical simulations. The final quantity 
is the quenched susceptibility at /3 = 5.93. 

We also expect, from the Maclaurin chiral expansion of Eqn. |20|, 

rVx = Co • (r M.) 2 - c 2 /c 3 ■ (f M.) 4 + c 3 /c 2 3 ■ (f M n ) 6 + O((f M n ) s ) (21) 

that c is related to the decay constant as before, c = (f / 7r ) 2 /4. We present the results 
of the fits in Table |. We find the UKQCD data to be well fitted by this form, but the 
asymptotic value is higher than the number we use for the quenched limit (in contrast to 

2 For an alternative motivation of this form, see fl23|. 
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earlier estimates of this value). Incorporating this number in the fit leads to a poorer x 2 , 
and a less robust fit. 



We turn now to fits based on the functional form in Eqn. [10]. We therefore use the 
ansatz 

where once again we expect C3 = ro 4 x qu and from the expansion 

r 4 x = CotfoikQ 2 - (^) 2 • (i) • (f M w ) 4 + O((r M„) 8 ) (23) 

we expect Co = {rofir) 2 /4. Note that in contrast to Eqn. [21], Eqn. |23] has no term that 
is cubic in m q and the rise will remain approximately quadratic for a greater range in 
(foMn) 2 . That this need be no bad thing is suggested by the relatively large range over 
which we could fit Eqn. [19|. Indeed, we see from the fits listed in Table || that this form 
fits our data quite well. 

Typical examples of the fits from Eqn. ^ and Eqn. [22] are shown in Figs. ^ and [5]. 
The similarity of the two functions is apparent. In Table |5] we use the fit parameters to 
construct the first three expansion coefficients in the Maclaurin series for the various fit 
functions, describing the chiral behaviour of x- The fits are consistent with one another. 

The fitted asymptote of the susceptibility at large m q is given by C3. We see from Table ^ 
that these are broadly consistent with the quenched value, and our large statistical errors 
do not currently allow us to resolve any 0(1/ TV 2 ) deviation from this. 

As an aside, we ask what happens if we cast aside some of our theoretical expectations 
and ask how strong is the evidence from our data that (a) the dependence is on M 2 rather 
than on some other power, and (b) the susceptibility really does go to zero as M n — > 0? 
To answer the first question we perform fits of the kind Eqn. but replacing (r M n ) 2 
by (f M 7r ) c . We find, using all seven values of x, that c = 1.32 (33) (10); a value broadly 
consistent with c = 2. The x 2 /d.o.f. is poorer, however, than for the fit with a power 
fixed to 2 (possible as there is one fewer d.o.f.) suggesting that the data does not warrant 
the use of such an extra parameter. As for the second question, we add a constant c 
to Eqn. ^ and find c = —0.056 (23) (25). Again this is consistent with our theoretical 
expectation; and again the x 2 /d.o.f. is worse. In both cases, however, the fits are not 
robust, with the fit parameters ill-constrained by our data. (See Table || for details of the 
above two fits.) 

Finally, we attempt to address the issue of discretisation effects. Our use of the non- 
perturbative c sw should have eliminated the leading 0(a) errors. Although (a/r ) 2 is 
small there is, of course, no guarantee that the coefficient of this correction might not 
be enhanced. We may begin to attempt to address this more quantitatively through a 
combined fit that includes the first order discretisation correction in tq. Our motivation 
here is not so much to give a continuum limit (our data will not really support reliably 
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such a long extrapolation in a 2 ) as to control the variations due to differing discretisation 
in our data over the relatively small range of lattice spacings in our study. For this reason 
we fit the simplest combined ansatz 

fix = co (r M n ) 2 + Cl (f M„) 4 + ^ (24) 

to the five most chiral of our data points. The resultant parameters are shown in Table [7|. 
The change in discretisation errors across out data set is clearly small, as we expected 
from the success of the previous fits. This justifies the use of a single chiral extrapolation 
over this limited range in fo- Whilst it does not, however, rule out deviations between 
the results of this and its equivalent in the continuum limit, it does give some indication 
that these deviations will be small. Clearly far greater accuracy in the measurements is 
needed to allow a confident extrapolation to a = 0. 

Given the consistency of our description of the small M x regime from our measurements, 
it is reasonable to use the values of Co to estimate the pion decay constant, f n . This is 
done in units of tq in Tables ^| and £|. We use the common chiral fit of Eqn. ^ over the 
largest acceptable range to provide us with our best estimate and its statistical error. We 
then use the fits with other functional forms to provide us with the systematic error. This 
produces an estimate 

r U = 0.262 ±0.015 IK (25) 

where the first error is statistical and the second is systematic. This is of course no 
more than our best estimate of the value of f n corresponding to our lattice spacing of 
a ~ 0.1 fm. This value will contain corresponding lattice spacing corrections and these 
must be estimated before making a serious comparison with the experimental value. We 
merely note that using f = 0.49 fm we obtain from Eqn. |25| the value 

U = 105 ±6 t\l MeV ( 26 ) 
which is reasonably close to the experimental value ~ 93 MeV. 



4 Comparison with other studies 



During the course of this work, there have appeared a number of other studies of the 



topological susceptibility in lattice QCD; in particular by the Pisa group |24|. 25(| , the CP- 
PACS collaboration p6|, |27|, the SESAM-T^L collaboration [28| and the Boulder group 

ESI. 



The most recent studies [27, 29] are consistent with our findings, but the earlier ones 
found no significant decrease of the susceptibility with decreasing quark (or pion) mass 
when everything was expressed in physical, rather than lattice, units. Indeed when our 
detailed results and analysis were first publicised || all the other studies then available 
23, (and indeed ]p8[| ) appeared to contradict our findings and it was therefore 
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necessary for us to provide some reason why this might be so. Although the situation is 
now different, the lessons are still useful and we will therefore briefly summarise the main 
point here. For more details we refer the reader to 0. 

All these other calculations differ from our study in having been performed at fixed /3. 
That implies that the lattice spacing a decreases as m q is decreased. In typical current 
calculations this variation in a is substantial. (See for example Fig. 4 of |3(J.) At the 
smallest values of m q the lattice spacing cannot be allowed to be too fine, because the total 
spatial volume must remain adequately large. This implies that at the larger values of m q 
the lattice spacing is quite coarse. Over such a range of lattice spacings, the topological 
susceptibility in the pure gauge theory typically shows a large variation (as, for example, 
in eqn.|T6|). Since for coarser a more instantons (those with p < 0(a)) are excluded, and 
more of those remaining are narrow in lattice units (with a correspondingly suppressed 
lattice topological charge) we expect that this variation is quite general, and not a special 
feature of the pure gauge theory with a Wilson action. In lattice QCD, we therefore 
expect two simultaneous effects in x as we decrease m q at fixed (3. First, because of the 
0(a 2 ) lattice corrections just discussed, x w iU (like x qu ) tend to increase. Second, it will 
tend to decrease because of the physical quark mass dependence. In the range of quark 
masses covered in current calculations this latter decrease is not very large (as we have 
seen in our work) and we suggest that the two effects may largely compensate each other 
so as to produce a susceptibility that shows very little variation with m q , in contrast to 
the ratio x/x qu which does. 

To illustrate this consider the fixed-/? calculation in The range of quark masses 
covered in that work corresponds to (r^M^) 2 decreasing from about 6.5 to about 3.0. 
Simultaneously a/ro decreases from about 0.437 to about 0.274. Over this range of 1/fo = 
a/r the pure gauge susceptibility increases by almost a factor of two, as we see using 
Eqn. [Tj| Clearly this is large enough to compensate for the expected variation of the 
susceptibility. 

As (3 is increased, the 0(a 2 ) variation with m q of the corresponding quenched suscep- 
tibility rpX qu will clearly diminish. Thus we suggested in || that if CP-PACS were to 
repeat their susceptibility calculations on their larger j3 ensembles, they would find an 
m q variation of TqX consistent with ours. This is what they have done in their most 
recent work [27|] with the result that we predicted. All this emphasises the utility of the 
UKQCD strategy of decoupling the variation of lattice corrections from the physical m q 
dependence, by performing calculations at fixed a rather than at fixed /3. 



5 Summary 

We have calculated the topological susceptibility in lattice QCD with two light quark 
flavours, using lattice field configurations in which the lattice spacing is approximately 
constant as the quark mass is varied. We find that there is clear evidence for the expected 
suppression of x with decreasing (sea) quark mass. 
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We have discussed this behaviour in the context of chiral and large N c expansions, and 
find good agreement with the functional forms expected there. We are not able to make 
a stronger statement about how close QCD is to its large iV c limit, owing to the relatively 
large statistical errors on our calculated values, particularly at larger quark masses. This 
situation should change in the near future and, together with the increasing availability of 
information on the large- N c behaviour of the pure gauge (quenched) theory 1TTJ], a more 
precise comparison will become possible. 

The consistent leading order chiral behaviour from our various fitting ansatze allows 
us to make an estimate for the pion decay constant, f % = 105 ± 6 jz^g MeV, for the 
lattice spacing of a ~ 0.1 fm. (Here the first error is statistical and the second has to 
do with the chiral extrapolation.) We use a lattice fermion action in which the leading 
0(a) discretisation errors have been removed. Since the more accurate (quenched) hadron 



masses show little residual lattice spacing dependence |^T], we might expect that this 
value of fn is close to its continuum limit. In any case, we note that it is in agreement 
with the experimental value, ~ 93 MeV. 
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label 




K 


Csw 


iVtraj. 


^0 


rh^/rhp 




ei 


5.20 


0.13565 


2.0171 


2800 


5.21 (5) 




1.386 (63) 


e 2 


5.20 


0.13550 


2.0171 


8300 


5.041 (40) 


0.5781^ 


1.480 (22) 


e 3 


5.25 


0.13520 


1.9603 


8200 


5.137 (49) 




1.978 (22) 


e 4 


5.20 


0.13500 


2.0171 


7800 


4.754 (40) 


0.7001$ 


1.925 (38) 


e 5 


5.26 


0.13450 


1.9497 


4100 


4.708 (52) 


0.783l| 


2.406 (27) 


e 6 


5.29 


0.13400 


1.9192 


3900 


4.754 (40) 


0.835t£ 


2.776 (32) 



Table 1: The ensembles studied, with measurements made every tenth HMC trajectory. 
Numbers in italics are preliminary estimates. 



ens. 


Tint(Q) 


(Q 2 ) 




ei 


32 (7) 


4.50 (72) 


0.0253 (41) 




48 (6) 


6.30 (96) 


0.0311 (48) 


e 3 




7.04 (66) 


0.0374 (37) 


e 4 


79 (5) 


11.48 (90) 


0.0447 (37) 


e 5 


89 (6) 


11.74 (1.00) 


0.0440 (42) 


e 6 


129 (19) 


17.13 (2.63) 


0.0668 (104) 


(3 = 5.93 




6.47 (24) 


0.049 (2) 



Table 2: The topological charge integrated autocorrelation time estimates (in units of 
HMC trajectories), the square of the charge and the susceptibility. The quenched result 
at /3 — 5.93 is for 16 4 fllT] , compared to 16 3 32 for the dynamical lattices. 



Fit 


iV fit 


c 


Cl 


X 2 /d.o.f. 


Q 




Eqn. 


2 


0.0141 (18) 




0.300 


0.584 


0.238 (16) 


Eqn. U 


3 


0.0106 (8) 




2.850 


0.058 


0.206 (8) 


Eqn. 


4 


0.0111 (7) 




2.311 


0.074 


0.211 (7) 


Eqn. H 


5 


0.0096 (5) 




5.289 


0.000 




Eqn. |IJ 


3 


0.0187 (36) (5) 


-0.0023 (11) (1) 


0.416 


0.519 


0.274 (27) 


Eqn. |19| 


4 


0.0184 (37) (5) 


-0.0020 (10) (1) 


1.405 


0.245 


0.272 (28) 


Eqn. p 


5 


0.0172 (19) (1) 


-0.0017 (4) (0) 


0.984 


0.399 


0.262 (15) 


Eqn. |19 


6 


0.0149 (16) (1) 


-0.0011 (3) (0) 


2.071 


0.082 


0.244 (14) 



Table 3: Fits to the iV fit most chiral points of (rlx)/M%. 
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Fit 
Eqn. |3 
Eqn. |g 
Eqn. |3 
Eqn. ^ 
Eqn. |3 



fit 



c 



3 
4 
5 
6 
7 



0.0262 
0.0216 
0.0257 
0.0220 
0.0445 



(142) (34) 
'29) 



(86 
(65 
(94 



24) 
20) 
21) 



0.0592 
0.0812 
0.0664 
0.0773 
0.0501 



c 3 

(252) (35) 
(379) (176) 
(156) (24) 
(202) (37) 
(20) (0) 



rVd.o.f. q 



0.553 
1.441 
1.123 
1.401 
2.134 



0.457 
0.237 
0.338 
0.231 
0.261 



0.324 (91) 
0.294 (64) 
0.320 (56) 
0.297 (46) 
0.422 (46) 



Eqn. 
Eqn. 
Eqn. 
Eqn. 
Eqn. 



3 
4 
5 
6 
7 



0.0210 
0.0190 
0.0205 
0.0184 
0.0219 



(69 
(52 
(44 
(37 
(30 



21) 
20) 
16) 
15) 
12) 



0.0456 
0.0564 
0.0517 
0.0575 
0.0497 



(115) (16) 
(158) (37) 
(73) (12) 
(95) (19) 
(19) (1) 



0.529 
1.408 
1.013 
1.457 
1.408 



0.467 
0.245 
0.386 
0.212 
0.300 



0.290 (51) 

0.276 (41) 

0.286 (34) 

0.272 (30) 

0.296 (22) 



Table 4: Fits to the N^ t most chiral points of (r^x). 



Fit 


iV fi t 


const. 


O((f M n ) 2 ) 




Eqn. m 


4 


0.0172 (20) 


-0.0017 (4) 




Eqn. |3 


5 


0.0220 (68) 


-0.0063 (43) 


0.0018 (20) 


Eqn. |3 


6 


0.0286 (39) 


-0.0130 (35) 


0.0059 (24) 


Eqn. |22| 


5 


0.0184 (40) 


-0.0024 (12) 





Eqn. 


6 


0.0170 (18) 


-0.0019 (4) 






Table 5: Chiral expansion terms of fitted functions. 



iVflt 


c 


c 3 


c 


c 


X 2 /d.o.f. 


Q 


7 


0.086 (32) (45) 


0.105 (23) (35) 


fixed to 2 


-0.056 (23) (35) 


1.570 


0.180 


7 


0.0160 (53) (26) 


0.0492 (20) (0) 


1.32 (33) (10) 


fixed to 


1.573 


0.178 



Table 6: Fits to the N& t most chiral values of (fjjx); for the power c of M|, and for the 
chiral intercept, c. These are as described in the text. 



N & t c ci 


c 2 


xVcl.o.f. 


Q 


5 0.0148 (75) -0.0014 (10) 


0.0055 (120) 


1.611 


0.200 



Table 7: Simultaneous chiral and lattice spacing correction fit to the iVg t most chiral 
values of (f^x), as described in the text. 
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Figure 1: The normalised variation of the topological charge with cooling, Rq(n c ) 
(Eqn. 0). 
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(3=5.20, k=0. 13550 




16000 18000 20000 22000 24000 
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Figure 2: A Monte Carlo time series of Q after ten cools. 
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p=5.20, K=0.13565 



P=5.26, K=0.13450 




Figure 3: Histograms of the topological charge for the most and second least chiral 
ensembles, together with the Gaussian given in Eqn. [15|. 
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Figure 4: The measured topological susceptibility, with quenched value at equivalent f . 
The radius of the dynamical plotting points is proportional to Tq . The fits, independent 



of the quenched points, are: (iii) Eqn. 20 and (iv) Eqn. g2| 
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Figure 5: The measured topological susceptibility. The radius of the dynamical plot- 
ting points is proportional to f$ . The fits, independent of the quenched points, are: 
(i) Eqn. [TJ, (ii) Eqn. ||, (iii) Eqn. || and (iv) Eqn. g| 
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